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We present a numerical study of the zero-temperature response of the Gaussian random-field Ising 
model (RFIM) to a slowly varying external field, allowing the system to be trapped in microscopic 
configurations that are not fully metastable. This modification of the standard single-spin-fiip 
dynamics results in an increase of dissipation (hysteresis) somewhat similar to that observed with 
a finite driving rate. We then study the distribution of avalanches along the hysteresis loop and 
perform a finite-size scaling analysis that shows good evidence that the critical exponents associated 
to the disorder-induced phase transition are not modified. 
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I. INTRODUCTION 

The random-field Ising model (RFIM) at zero tempera- 
ture has been proposed as a prototype for a broad class of 
disordered systems (random magnets, glasses, plastic and 
ferro-elastic materials, fluids in porous media ...) which 
exhibit an intermittent, avalanche- like, and hysteretic re- 
sponse to a smoothly varying applied force The RFIM 
has also been used in several socio-economics contexts to 
simulate collective effects induced by imitation and so- 
cial pressure [1]. A remarkable prediction of the model 
is the existence of a non-equilibrium critical point (for a 
certain amount of disorder) which separates two differ- 
ent regimes of avalanches. In the strong-disorder regime, 
all avalanches are of microscopic size and the saturation 
hysteresis loop (e.g. magnetization m versus magnetic 
field H) is macroscopically smooth. At low disorder, a 
macroscopic avalanche occurs at a certain field, which re- 
sults in a jump discontinuity in the magnetization. At 
criticality, avalanche sizes and durations follow power-law 
distributions. 

In the original version of the model |3|, spins obey 
a standard single-spin-flip (Glauber) relaxation dynam- 
ics at T = and align with their local effective field. 
As the applied field is slowly increased or decreased, a 
spin may become unstable and then trigger an avalanche 
that propagates until another metastable state is found. 
The field is held fixed during the propagation of the 
avalanche, which corresponds to the so-called "adiabatic" 
limit. This amounts to assuming that the time for lo- 
cal equilibration (the duration of an avalanche) is much 
smaller than the rate of change of the driving field. More- 
over, by using a deterministic zero-temperature dynam- 
ics, one assumes that no thermally activated escape takes 
place on the observational time scale ("athermal" limit). 
These two assumptions are reasonable in many physical 
situations but they are never fully satisfied. It is thus in- 
teresting to test the robustness of the predicted scenario 



(in particular the existence of a disorder-induced phase 
transition and its universality) with respect to a slight 
violation of these conditions. 

In a recent workQ, the influence of partial equilibra- 
tion processes was mimicked by changing the dynamics 
and allowing two neighboring spins to flip cooperatively. 
As expected, this change resulted in a reduction of hys- 
teresis (as the set of two-spin-flip stable states is con- 
tained in the set of one-spin-flip stable states) and an 
enhancement of large-scale collective effectsQ. But, re- 
markably, the critical behavior (characterized by vari- 
ous exponents and finite-size scaling functions) remained 
identical. In the present work, we want to go in the oppo- 
site direction by enlarging the set of visited microscopic 
states so to increase hysteresis and drive the system fur- 
ther from equilibrium. This is done by allowing a small 
fraction of the spins to point in the opposite direction 
to their local field. However, in order to keep things as 
simple as possible (and, in particular, to keep the simpli- 
fying separation of motion between adiabatic driving and 
avalanche propagation) , we will still use a single-spin-flip 
dynamics and start an avalanche when some threshold in 
energy is reached. As a consequence, the saturation loop 
will not be an "extremal" path in the field-magnetization 
plane^G^, 't'I and the property of return-point memory}^ 
will be violated. The main question that we want to 
address is: will this change the universal properties of 
the critical behavior? Note that this modification of the 
dynamical rule may be viewed as a crude way of simu- 
lating the effect of a finite driving rate which does not 
give enough time to the system to relax, even locally 
[1, [^, [13, [m, [13 • But this is of course a caricature of 
what happens in the real world and we therefore do not 
pretend to propose here a theory of dynamic hysteresis, 
a topic that has been (and still is) intensively investi- 
gated in the literature fl3j. In a socio-economic context, 
one could also interpret this model as simulating an ef- 
fect of "inertia" that prevents the individual agents to 
make their decision (e.g., buy or sell) when the incentive 
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reaches the threshold value. 

The paper is organized as follows. In Section|TT], we de- 
scribe the model and the new dynamics, and discuss some 
of its properties. In Section [1111 we study the hysteresis 
loops, in particular the change in coercivity. Avalanche 
distributions are analyzed in Section IIVI and the criti- 
cal behavior in Section |Vl Summary and conclusion are 
given in Section IVTl 

II. MODEL AND DYNAMICS 

We consider the RFIM on a three-dimensional cubic 
lattice of linear size L with periodic boundary conditions. 
On each site i there is a spin variable Si = ±1. The 
energy of the system of = L'^ spins is described by the 
Hamiltonian 

n^^JY,s,s,-Y,h^s^~HY,s^, (1) 

<ij> i i 

where J > is a ferromagnetic exchange interaction con- 
stant and the first sum runs over nearest-neighbor pairs. 
Without loss of generality, we will take J = 1. The 
set of random fields {hi} is drawn independently from 
a Gaussian probability distribution with zero mean and 
variance a^, and _ff is a uniform external magnetic field 
which couples to the overall magnetization M = Si 
(m = M/N). 

The energy change associated to the reversal of spin i 
then reads 

AH, = An{S^ ^ -S,) = 2S^F„ (2) 

where 

j/i 

is the local effective field acting on spin i, and j is a 
nearest- neighbor of i. 

At T = the standard single-spin-fiip dynamics con- 
sists in fiipping a spin if this lowers its energy. A con- 
figuration is thus single-spin-flip stable if every spin is 
aligned with its local field, i.e., 

S^^sign{F,) Vi. (4) 

Each of these metastable states has a certain range of 
stability 7J,„j„ < -ff < H„iax- 

As the field H is slowly changed, a spin flips (either 
upwards or downwards) when its local field changes sign, 
which may induce a whole avalanche of other spin flips. 
In the adiabatic limit H is held constant during the prop- 
agation of the avalanche. The avalanche stops when a 
new metastable state is reached. A nice feature of this 
dynamics is that it is "abelian" : because of the so-called 
no-passing " ruled, the metastable state reached at the 
end of an avalanche is independent of the order in which 
unstable spins have been reversed. In consequence, these 



spins can be flipped either sequentially or in parallel (this 
latter choice having the advantage that one can define the 
duration of an avalanche). 

We now modify this dynamical rule by allowing a cer- 
tain number of spins to be unstable (Eq. (j4]) is then vi- 
olated). This can be done in various ways, for instance 
by imposing that the fraction of unstable spins cannot 
exceed a certain value. This, however, would prevent 
the system to reach saturation. We therefore prefer to 
compute the extra amount of energy associated to the 
unstable spins and impose that an avalanche starts when 
this contribution exceeds some fixed threshold e. Specif- 
ically, the system may visit spin configurations that we 
call "e-stable" and that satisfy 

J2 AH, > -^e, (5) 

Si unstable 

where the sum runs over all unstable spins, i.e. spins for 
which ATii < (note that e is an intensive quantity). 
There are now two conditions for an avalanche to start: 
(i) there must be unstable spins (ii) the sum of the ex- 
tra contributions to the energy due to the unstable spins 
must exceed the threshold (in absolute value). When in- 
creasing (resp. decreasing) the field, only negative (resp. 
positive) spins contribute to this energy. One of course 
recovers the usual dynamics for e = 0. 

It is quite obvious that this modification of the dynam- 
ics spoils the no-passing rule and the abelian property. 
Therefore, since several spins may be unstable when the 
energy threshold is reached, one needs to specify the or- 
der in which these spins are flipped. The natural choice 
that we adopt is to flip the most unstable spin first, that is 
the one that corresponds to the most negative AHi. We 
then update the local fields while keeping H constant, 
search again for the most unstable spin, fiip it, etc., until 
Eq. [5] is satisfied [U]. This "greedy" algorithm essentially 
amounts to performing a steepest descent path in energy, 
which defines a deterministic sequence of unstable states 
inside an avalanche (defined as usual as the collection of 
spins which flip at the same field) . It is then not difficult 
to see that the order in which the spins are fiipped along 
the hysteresis loop does not depend on e, so that the 
sequence of visited states is invariant. Changing e only 
changes the values of the external field at which the states 
are reached or left. More generally, in any monotonous 
field history from a state with magnetization M to a state 
with magnetization M' , the sequence of visited spin con- 
figurations does not depend on e, but the corresponding 
fields may differ. As a result, avalanches may split or 
merge when changing e so that their number and size 
also change. This is illustrated in Fig. [T] where we com- 
pare the hysteresis loops obtained for e = and e = 0.03 
in a small size system. One can see in Figs. [Ijb) and 
(c) that the spin configurations at a given magnetization 
M are identical although the hysteretic trajectories are 
different. Fig. [Hb) also shows that a state which is not 
single-spin-flip stable for e = (as it is located in the 
middle of an avalanche) is "e-stable" for e = 0.03. In 
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FIG. 2: (Color online) Test of the return-point memory prop- 
erty in a system of linear size L = 12 for cr = 3.5 and e = 0.05. 
The field is increased up to Hq, decreased from Hq to H\ and 
then increased again up to saturation. The return trajectory 
crosses the ascending branch of the saturation hysteresis loop 
several times, indicating that e-states exist outside the loop 
and that the RPM property is violated. 



FIG. 1: (Color online) Comparison of the hysteresis loops ob- 
tained with e = (solid line) and e = 0.03 (dashed line) in 
a system of linear size L = 12. In (b) and (c) the upper and 
lower parts of the ascending branches are blown up, show- 
ing that some avalanches merge or split when e is changed. 
For a given overall magnetization M, the reached configura- 
tions are identical (as illustrated by two-dimensional slices of 
the system where negative spins are drawn in black), but the 
corresponding fields are different. 

Fig. [He), the opposite situation is observed. 

The fact that the modified dynamics does not satisfy 
the no-passing rule has two consequences. Firstly, there 
can exist e-states outside the hysteresis loopf?]. Secondly, 
the property of return-point memory (RPM) is not satis- 
fied, as shown in Fig. [2] However, the violation of RPM 
is small for the values of e considered in this work and 
it seems that this property is better and better verified 
as the system size increases. One can also see in Fig. [5] 
that the first-order reversal curves have a linear portion: 
when reversing the field (for instance from 7?o to Hi)^ 
the negative spins that were unstable at i?o become sta- 
ble again before any positive spin becomes unstable and 
Eq. [5]is only violated when H < 7?2- 

As a final remark in this Section, we compare the pro- 
posed dynamics with the algorithms that have been used 
previously to study the T=0 RFIM with a finite driving 
rate s, IQ, [H 

[T^ . The most common strategy consists 
in increasing the field in finite steps AH, thus merging 
all the avalanches occurring within that field window into 
a larger avalanche [1, [ll|, [l3]. The resulting magnetiza- 
tion loops then share a series of common points with 
those obtained with the adiabatic driving. On the other 



hand, in our case avalanches not only merge but also 
split and the resulting loops differs everywhere from the 
adiabatic (e = 0) loops. A second strategy [l3, [HI con- 
sists in performing an exact simulation of the continuous 
M{t) signal by using a finite driving rate and defining a 
time interval associated to the shell of spins that relax in 
parallel. But one then needs to fix a threshold to define 
the avalanches and this has a strong influence on their 
size[3|. This problem does not occur in the modified 
dynamics that we used here since there is still a com- 
plete time-scale separation between the field driving and 
the propagation of the avalanches. 



III. HYSTERESIS LOOPS 

We first consider the infiuence of e on the shape of 
the hysteresis loops for different values of the disorder 
strength a. Fig. [3] shows the results obtained by aver- 
aging over different disorder realizations for a — 3 and 
a = 1.5. As expected, the main effect of e is to bring the 
system further away from equilibrium and to increase 
hysteresis. For a — 3 for instance, the loop area, which 
represents the energy loss, increases by ~ 70% when in- 
creasing e from to 0.03. This is already an important 
variation and in the following we shall restrict our study 
to the range < e < 0.03. 

Fig. [3] also shows that there are still two different 
regimes when e ^ and that the shape of the loops 
changes from smooth to rectangular as a decreases. In 
particular, at low disorder, there is a single avalanche 
that spans the whole system at a certain value of the ex- 
ternal field (note that the spanning avalanche in a finite 
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FIG. 3: (Color online) Average hysteresis loop for a — 2.8 
(a) and a = 1.5 (b). The curves correspond respectively to 
e = (inner loop), 0.0005,0.004,0.008,0.01,0.02, and 0.03 
(outer loop). Data result from an average over 5000 disorder 
realizations in a system of size L — 24. 

system occurs with no collective precursor for e > 0.03, 
which may change its nature. This is also a reason to 
restrict our study to smaller values of e). 

To further quantify the influence of e on the hysteresis, 
we show in Fig. S] the variation of the coercive field Hcoer 
for different values of the disorder. We find that the data 
are quite accurately fitted by the equation 

Hcoerie) = HcoeriO) + Ce^ , (6) 

with P « 0.5 in the large-disorder regime and /3 w 0.45 
at low disorder [l^. The same dependence is found for 
the variations of the loop area with e. So far, we have 
no convincing theoretical explanation for this behavior. 
On the other hand, it appears that the same exponent 
P « 0.45 in the low disorder regime has been observed 
in a simulation study of the RFIM under a linear driving 
ratefoll (in that study, however, there is no clear scaling 
at large disorder). In two dimensions (l^ . when varying 
the field in small steps as in Ref. |Tl|, simulations show 
a crossover from a square-root to a linear dependence 
of Hcoer with the rate as the disorder is increased, in 
agreement with the behavior observed in ferromagnetic 
thin films. 



IV. AVALANCHE SIZE DISTRIBUTION 

The RFIM with the standard (e = 0) dynamics dis- 
plays a power-law distribution of avalanche sizes at a 
critical disorder cr = CTc[3]. Depending on the method 
used to extrapolate the numerical results to the thermo- 
dynamic limit, the value of CTc varies from 2.16 [17] to 
2.21 [H, U^. In this section we study the behavior of the 
avalanche size distribution for e > 0. 

Fig. [5] shows the avalanche size distribution 
D{s; a, e, L) obtained in a system of size L = 30 
for e = 0.008 and various disorders a when sweeping 



FIG. 4: Average coercive field as a function of e for different 
values of a. The data result from an average over 500 disorder 
realizations in a system of size L = 24 and are fitted according 
to Eq. (O. 
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FIG. 5: (Color online) Avalanche size distributions in a sys- 
tem of size L = 30 for e — 0.008 and different values of a. 
The curves are sorted from top to bottom in the order indi- 
cated in the legend. The statistics has been performed over 
1500 disorder realizations. Data for a 7^ 2.3 have been shifted 
vertically. 

through a half-loop. The same behavior as for e = is 
observed: for large a the distribution is exponentially 
damped whereas for small a there is a peak at large sizes 
due to avalanches with characteristic size ^ L^. These 
spanning avalanches are responsible for the macroscopic 
discontinuity in the hysteresis loop. Between these two 
regimes, there is a value of a for which the distribution is 
very well approximated by a power-law up to the trivial 
cutoff size Smax ^ (and apart from some corrections 
at very small s). This "critical" value of a changes with 
e but the slope of the power-law region appears to be 
invariant, as shown in Fig. [SI. We view this result as a 
first indication that there exists a critical disorder crc(e) 
in the thermodynamic limit and that the power-law 
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exponent t + cr[35 that characterizes the avalanche size 
distribution at criticaHtyQ does not change with e (note 
that the numerical value obtained in a finite system may 
differ from the actual value ~ 2 in the thermodynamic 
limit, as estimated in Ref. [^. In the next Section we 
present a finite-size scaling analysis that will corroborate 
this statement. 

Note that this result contrasts with the one that has 
been obtained in previous studies using a finite driving 
rate [l^, [IH ■ Indeed, when the only effect of the finite 
driving rate is to merge small avalanches into larger ones, 
the power-law exponent decreases. It is in fact unclear 
if such a decrease is due to the actual out-of-equilibrium 
behavior or is induced by the approximate treatment of 
the avalanche merging phenomenon and/or by the def- 
inition of the threshold that allows to discriminate the 
avalanches. 
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FIG. 7: (Color online) Average number of ID- (a) and 2D- 
(b) spanning avalanches as a function of a for e = 0.008 and 
different system sizes. Typical error bars for the largest sizes 
are shown. 



FIG. 6: (Color online) Avalanche size distributions in a sys- 
tem of size L = 30 for different values of e. The curves are 
sorted from top to bottom in the order indicated in the legend. 
In each case, the disorder a is the one for which the distribu- 
tion is closest to a power-law. The slope of the dashed lines 
that describe the power-law region is —1.8. The statistics 
has been performed over 1500 disorder realizations. Data for 
e 7^ 0.02 have been shifted vertically. 



V. CRITICAL PROPERTIES 

The analysis of the spanning avalanches has proven 
to be a successful way to determine the values of sev- 
eral critical exponents in the RFIM with the standard 
metastable dynamics [l^, [l^. Indeed, the statistics of 
spanning avalanches in finite systems contains informa- 
tion about the percolating fractal avalanches in the ther- 
modynamic limit which are the signature of criticality. In 
this Section we use a finite-size scaling method to study 
the effect of e on the critical behavior of the model. 

With the modified dynamics (e > 0) an avalanche 
still involves a connected set of spins and the spanning 



avalanches can be defined as usual : avalanches that span 
the whole system in 1,2, or 3 spatial dimensions are re- 
ferred to as ID-, 2D-, and SD-spanning avalanches, re- 
spectively. We focus our analysis on the average num- 
ber of ID- and 2D-spanning avalanches occurring along 
the lower branch of the hysteresis loop, which we call 
7Vi(cr, e, L) and N2[<J, e, L), respectively. We discard from 
our analysis the 3D-spanning avalanches, which con- 
tain information not only about the critical percolating 
avalanches, but also about the compact infinite avalanche 
that gives rise to the first-order discontinuity in the low- 
disorder regime [ill . 

Figure [7] illustrates the behavior of Ni{a,e,L) and 
N2{a,e,L) as a function of cr for e = 0.008 and differ- 
ent system sizes (the data correspond to averages over 
typically 10^ — 10'' disorder realizations). One can see 
that both functions exhibit a peak whose height increases 
with L. Moreover the peak position shifts and its width 
reduces. These features are clear signatures of the ex- 
istence of an infinite number of percolating avalanches 
at a certain critical disorder (Tc(e) in the thermodynamic 
limit. Ni{a^ e, L) and N2{u^e,L) are thus expected to 
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FIG. 8: (Color online) Height and position of the peak in 
A''i (cr, £, L) and N2 (cr, e, L) for different values of e and differ- 
ent system sizes. In (a), the height is plotted vs. L in a log-log 
plot and in (b) the position is plotted vs. L"^''" with v — 1.2. 
The data are fitted according to Eqs.lO with A — —0.1. The 
dotted line in (a) has a slope 0.1 and is included as a reference. 



have the scalmg formp^Tl iTsj 

iVi(a,e,L) = L^iVl(wil/^e), (7) 
N2{a,e,L) = MuL^'" .e) , (8) 

where A'^i and N2 are finite-size scaling functions, u{a) 
is some analytical function of the distance to the crit- 
ical disorder, and 6 and v are critical exponents that 
characterize the divergence of the number of spanning 
avalanches and of the correlation length, respectively. Al- 
though the function w(cr) can be approximated to first- 
order as (cr — (Tc(e))/fc(e), previous studies pjj] suggest 
that higher order terms are necessary to produce good 
scaling collapses. Therefore, as in Ref. [H, we shall use 



g - crc(e) 
(Tc(e) 



A{e) 



crc(e) 



(Tc(e) 



(9) 



where A{e) is a nonuniversal parameter that may depend 
on e. For e = 0, the best choice was found to be ^ = 
—0.2. This value did not change when replacing the 1- 
spin-flip by the 2-spin-fiip dynamics Q- 




FIG. 9: (Color online) Scaling plot of the number of ID- and 
2D-spanning avalanches for e = 0.008 and different system 
sizes. The upper (resp. lower) curve corresponds to the ID- 
(resp. 2D-) spanning avalanches using tJc = 1.97, A = —0.1, 
e = 0.1, and V = 1.2. 



Equations ([7]) and ([H]) may be first tested by plotting 
the height of the peaks in A^i(cr, e, L) and N2{(y,e,L) as 
a function of L in a log-log scale. As shown in Fig.[8ja), 
the behavior is then linear and the slope is compatible 
with the value 9 = 0.1 obtained for e = ofTsj. 

A second test consists in plotting the position cr„iax of 
the peaks as a function of L^^/'^ , as shown in Fig. Mjo)- 
From Eqs. ([7]) and ([8]), this position should be determined 
by the condition uL^^/'^ — constant, i.e.. 



,(e) - (Jc{e) 
crc(e) 



A 



,(e) -ac{e) 
CTc(e) 



(10) 

For yl = this equation predicts a linear behavior of 
(^max as a function of i^^/". As can be seen in Fig.[5I^b), 
the actual behavior is indeed almost linear when using 
the value v = 1.2 obtained for e = o[2| and the data for 
the ID- and 2D-spanning avalanches reasonably extrap- 
olate towards the same value in the thermodynamic limit 
L 00. This method, however, cannot be used to ex- 
tract accurate values of A or (e) and the fits shown in 
Fig. Eljb) are based on the results of the finite-size scaling 
analysis that we now discuss. 

Indeed, the best way to estimate all universal and 
nonuniversal parameters is to search for a good collapse 
of the curves Ni{a, e, L) and N2{(7, e, L) corresponding to 
different sizes L. This is done by plotting Ni{a, e, L)L~^ 
as a function of the scaling variable uL^/^ . As an ex- 
ample, we show in Fig. [9] the best collapse obtained for 
e = 0.008 using the values CTc = 1.97, A = -0.1, 9 = 0.1, 
and v — 1.2. This procedure can be done independently 
for each set of data corresponding to different e. Table 
U shows the parameters that produce the best collapses 
of the curves for 0.0005 < e < 0.03(2^ For comparison 
we also include the results obtained for e = ofl^. In all 
cases, the values v = 1.2, 9 = 0.1, and A = —0.1 are the 
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TABLE I: Universal and nonuniversal parameters that yield 
the best finite-size scaling collapses for different values of e. 



^ 0.6 
0.4 



0.2 ^ 



L 


= 12 


■ 


L 


= 16 


• 


L 


= 20 


A 


L 


= 24 


T 


L 


= 30 


♦ 




10 12 



best estimates. The only parameter showing a clear de- 
pendence with e is (Tc- Notice that A, which in principle 
is a nonuniversal parameter, takes the same value —0.1 
for all e > 0. For e = the value A = —0.2 produced a 
better collapse [l^ but the difference is not very signifi- 
cant since the collapse in Fig. [9] is still rather good when 
using this value (alternatively, one can also use A = —0.1 
for e = 0). 

The most important feature in Table [J is that the two 
critical exponents v and 9 do not change with e and are 
the same as for e = 0. This suggests that the critical 
behavior is described by the same universality class. We 
then try to collapse all the data for different e and L on 
the same plot by introducing a nonuniversal scale factor 
B that is e-dependent, i.e., by assuming that 

iVi(a,e,L) = L''Ni{B{t)uL^''') (11) 
N2{a,e,L) = L' N2{B{e)uL^''') . (12) 

As can be seen in Fig. [TUl a very good collapse of the 
whole set of curves in the range 0.0005 < e < 0.03 is 
obtained with the values of B indicated in Table |T] (this 
scale factor is arbitrarily set equal to 1 for e = 0.008). 

A cross-check of the consistency of the scaling collapses 
can be done by fitting the data in Fig.[5I^b) using Eg. [TUl 
with V = 1.2 and A = —0.1. The extrapolated values of 
CTc for L ^ oo are fully compatible with those reported 
in Table n 

Our data are thus consistent with the fact that the 
disorder-induced critical point found for e = transforms 
into a critical line when e > and the system is allowed 
to visit weakly unstable states. The whole critical line 
appears to be described by the same exponents and by 
the same scaling function for e > 0. On the other hand, 
it seems that the scaling function differs from the one 
corresponding to e = 0, as shown in Fig. [Tl] on a linear- 
log scale (in the figure, B is chosen so as to produce 
the best collapse on the right-hand side of the peak; for 
e = 0, this yields B = 1.35 and A = —0.1 but the picture 
essentially does not change with A = —0.2). 

The fact that the exponents v and 9 are the same but 
the scaling functions are different for e = and e > 
may be surprising at first sight. On the one hand, this 



FIG. 10: (Color online) Scaling plot of the number 
of ID- (upper curves) and 2D- (lower curves) spanning 
avalanches according to Eqs. (|lip and (|12|) for different sys- 
tem sizes L (indicated by different symbols) and for e = 
0.001, 0.002, 0.006, 0.008, 0.01, 0.02, 0.03 (not indicated). 




FIG. 11: (Color online) Comparison of the scaling functions 
iVi (above) and 7^2 (below) for e = (empty symbols) and 
e > (filled symbols) on a linear-log scale. Different symbols 
indicate the system size and different colors the value of e. N2 
has been shifted two decades downwards for clarity. 



may simply indicate that higher-order terms in the scal- 
ing variable u are needed (let us recall again that the 
simplest choice u — {a — o'c)/fc does not produce good 
scaling collapses and that it was proposed already in 
Ref. 3 to use u — [a — ac)/cr, which amounts to keep- 
ing an infinite number of terms in an expansion in pow- 
ers of ((T — (Tc)/(Tc). It is clear that simulations with 
much larger system sizes would be required to fully set- 
tle this issue[2l[. On the other hand, there may be in- 
deed an essential physical difference in the properties of 
the percolating clusters at the critical point. This may 
not change the fractal dimension (as this would be re- 
flected in a change of the critical exponents) but only the 
way the finite size affects the number and size of these 
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clusters. It can be seen for instance in Fig. [TT] that the 
number of ID- and 2D-spanning avalanches diverges like 
^ L^-^ at the critical point u — Q but that the prefac- 
tor of this divergence dramatically decreases as soon as 
e becomes slightly positive. Finite-size scaling functions 
are known to depend sensitively on boundary conditions 
and, in the present case, they may account for the fact 
that the percolating clusters are growing in a different 
environment when e > 0. Recall that we are dealing here 
with a nonequilibrium phenomenon and that, in the orig- 
inal model[l4|, a finite fraction of the system has already 
transformed when the critical point is reached. This is 
reflected in the value of the critical magnetization Mc- It 
is not easy to estimate the actual value of this quantity 
in the thermodynamic limit, but preliminary calculations 
suggest that Mc significantly decreases when e > (even 
for e = 0.0005), showing that the fraction of spins that 
flip when driving the system from H = — oo to the critical 
held i?c(e) becomes very small. 

VI. SUMMARY AND CONCLUSION 

We have studied the zero-temperature random-field 
Ising model, with a Gaussian distribution of the random 
fields, using a modified single-spin-fiip dynamics that al- 
lows the system to be trapped in weakly unstable states 



when driven quasi-statically by an external field. The 
new dynamics, however, does not modify the sequence of 
states that are visited during a monotonous field history 
and preserves the intermittent, avalanche-like character 
of the response to the driving field. The violation of the 
standard local stability condition is controlled by a single 
parameter e whose effect is somewhat similar to that of a 
finite-driving rate, moving the system away from equilib- 
rium and resulting in a similar increase of the width of the 
saturation hysteresis loop, as measured by the coercive 
field. Avalanches are modified but two distinct regimes 
of avalanches and two different loop shapes as a function 
of disorder are still present. As in the original model[^, 
the transition between the two regimes corresponds to a 
critical point where avalanches of all sizes are observed. 
The critical exponents that have been extracted from a 
finite-size scaling analysis of the number of the spanning 
avalanches appear to be independent of e, suggesting that 
the condition of strict metastability (as well as the no- 
passing rule) may be irrelevant for the critical behavior. 
In our opinion, this significantly enlarges the domain of 
validity of the original model. 
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